Nearly Periodic Maps and Geometric Integration of Noncanonical Hamiltonian Systems

M. Kruskal showed that each continuous-time nearly periodic dynamical system admits a formal U(1)-symmetry, generated by the so-called roto-rate. When the nearly periodic system is also Hamiltonian, Noether’s theorem implies the existence of a corresponding adiabatic invariant. We develop a discrete-time analog of Kruskal’s theory. Nearly periodic maps are defined as parameter-dependent diffeomorphisms that limit to rotations along a U(1)-action. When the limiting rotation is non-resonant, these maps admit formal U(1)-symmetries to all orders in perturbation theory. For Hamiltonian nearly periodic maps on exact presymplectic manifolds, we prove that the formal U(1)-symmetry gives rise to a discrete-time adiabatic invariant using a discrete-time extension of Noether’s theorem. When the unperturbed U(1)-orbits are contractible, we also find a discrete-time adiabatic invariant for mappings that are merely presymplectic, rather than Hamiltonian. As an application of the theory, we use it to develop a novel technique for geometric integration of non-canonical Hamiltonian systems on exact symplectic manifolds.


Introduction
A continuous-time dynamical system with vector parameter γ is nearly periodic if all of its trajectories are periodic with nowhere-vanishing angular frequency in the limit γ → 0. Examples from physics include charged particle dynamics in a strong magnetic field, the weakly relativistic Dirac equation, and any mechanical system subject to a high-frequency, time-periodic force. In the broader context of multi-scale dynamical systems, nearly periodic systems play a special role because they display perhaps the simplest possible non-dissipative short-timescale dynamics. They therefore provide a useful proving ground for analytical and numerical methods aimed at more complex multi-scale models.
In a seminal paper (Kruskal 1962), Kruskal deduced the basic asymptotic properties of continuous-time nearly periodic systems. In general, each such system admits a formal U (1)-Lie symmetry whose infinitesimal generator R γ is known as the rotorate. In the Hamiltonian setting, existence of the roto-rate implies existence of an all-orders adiabatic invariant μ γ by way of Noether's theorem. General expressions for μ γ may be found in Burby and Squire (2020). Recently (Burby and Hirvijoki 2021), we extended Kruskal's analysis by proving that the (formal) set of fixed points for the roto-rate is an elliptic almost invariant slow manifold. Moreover, in the Hamiltonian case, we demonstrated that normal stability of the slow manifold is mediated by Kruskal's adiabatic invariant. The purpose of this article is to introduce discrete-time analogs of continuoustime nearly periodic systems that we call nearly periodic maps. These objects can be motivated as follows. A nearly periodic system characteristically displays limiting short-timescale dynamics that ergodically covers circles in phase space. This ergodicity is ultimately what gives rise to Kruskal's roto-rate and, in the presence of Hamiltonian structure, adiabatic invariance. It is therefore sensible to regard parameter-dependent maps whose limiting iterations ergodically cover circles as discrete-time analogs of nearly periodic systems. Ergodicity requires that the rotation angle associated with each circle be an irrational multiple of 2π . In principle, these rotation angles could vary from circle to circle, but smoothness removes this freedom and imposes a common rotation angle across circles. Nearly periodic maps are defined by limiting iterations that rotate a family of circles foliating phase space by a common rotation angle. Such a map is resonant or non-resonant when the rotation angle is a rational or irrational multiple of 2π , respectively. The preceding remarks suggest that non-resonant nearly periodic maps should share important features with continuous-time nearly periodic systems.
We will show that non-resonant nearly periodic maps always admit formal U (1)symmetries by modifying Kruskal's construction of a normal form for the roto-rate. Thus, non-resonant nearly periodic maps formally reduce to mappings on the space of U (1)-orbits, corresponding to elimination of a single dimension in phase space. In the Hamiltonian setting, we will establish a discrete-time analog of Noether's theorem that will allow us to construct all-orders adiabatic invariants for non-resonant nearly periodic maps. In contrast to the continuous-time case, there may be topological obstructions to the Noether theorem-based construction. Nevertheless, assuming (a) existence of a fixed point for formal U (1)-symmetry, or (b) existence of a timedependent Hamiltonian suspension for the nearly periodic map ( a time-dependent Hamiltonian flow that interpolates between the identity and the nearly periodic map), these topological obstructions disappear. When an adiabatic invariant does exist, the phase-space dimension is formally reduced by two instead of one. On the slow manifold, corresponding to vanishing of the adiabatic invariant, the reduction in dimensionality may be even more dramatic.
We anticipate that non-resonant nearly periodic maps will have important applications to numerical integration of nearly periodic systems. While development of integrators for such systems is straightforward when the numerical timestep h resolves the short-timescale dynamics, considerably more care is required when "stepping over" the period of limiting oscillations. One approach would be to design an integrator on the unreduced space that is constrained to be a non-resonant nearly periodic map. Although such an integrator would not accurately resolve the phase of short-scale oscillations when taking large timesteps, it would automatically possess an all-orders reduction to the space of U (1)-orbits. By designing the reduced map to discretize the continuoustime reduced dynamics, the slow component of the continuous-time dynamics could be accurately resolved without directly simulating the reduced dynamical variables. This opens the door to a type of asymptotic-preserving integrator capable of seamlessly transitioning between large-and small-timestep regimes, generalizing those proposed in Ricketson and Chacón (2020), Xiao and Qin (2021) for magnetized charged particle dynamics. Moreover, in the Hamiltonian case, the integrator would automatically enjoy an all-orders adiabatic invariant close to the continuous-time invariant. Such a capability would complement previous results on short-timestep adiabatic invariants for variational integrators (Hairer and Lubich 2020). We provide a proof-of-principle demonstration of these ideas in Sect. 5.1 Aside from serving as integrators for nearly periodic systems, nearly periodic maps may also be used as tools for structure-preserving simulation of general Hamiltonian systems on exact symplectic manifolds. (See Abraham and Marsden 2008;Marsden and Ratiu 1999 for the foundations of Hamiltonian mechanics on symplectic manifolds.) The basic idea is to first embed the original Hamiltonian system as an approximate invariant manifold inside of a larger nearly periodic Hamiltonian system, as discussed in Burby and Hirvijoki (2021). Then, it is possible to construct a symplectic nearly periodic map that integrates the larger system while preserving the approximate invariant manifold. Discrete-time adiabatic invariance ensures that the approximate invariant manifold enjoys long-term normal stability, which is tantamount to the integrator providing a persistent approximation of the original system's dynamics. We describe and analyze this construction in Sect. 4.2. In Sect. 5.2, we apply the general theory to the non-canonical Hamiltonian dynamics of a charged particle's guiding center (Northrop 1963;Littlejohn , 1983) in a magnetic field of the form B = B(x, y) e z (Littlejohn 1979).
The remainder of this article is organized as follows. After providing a brief non-technical overview, we review Kruskal's theory of nearly periodic systems using modern terminology in Sect. 3. Then, we develop the general theory of nearly periodic maps in Sect. 4, including their special properties in the symplectic case, and their ability to serve as geometric integrators for Hamiltonian systems on exact symplectic manifolds. Wherever possible, proofs of general properties of nearly periodic maps parallel Kruskal's arguments from the continuous-time setting. Section 5 contains a pair of interesting applications of nearly periodic map technology. Finally, Sect. 6 provides additional review and context for this work.

Notational Conventions
In this article, smooth shall always mean C ∞ , and will always denote a vector space. We reserve the symbol M for a smooth manifold equipped with a smooth auxiliary Riemannian metric g. We say f γ : M 1 → M 2 , γ ∈ , is a smooth γ -dependent mapping between manifolds M 1 , M 2 when the mapping is an element of the tensor algebra T m (M) at m for each m ∈ M and γ ∈ , and (b) T γ is a smooth γ -dependent mapping between the manifolds M and The symbol X γ will always denote a smooth γ -dependent vector field on M. If T γ is a smooth γ -dependent section of either T M ⊗ T M or T * M ⊗ T * M, then T γ is the corresponding smooth γ -dependent bundle map T * M → T M : α → ι α T γ , or T M → T * M : X → ι X T γ , respectively. Note that if is a symplectic form on M with associated Poisson bivector J then −1 = − J . Finally, we introduce the following definition to address presymplectic forms that depend on a parameter.

Definition 1 Let
γ be a finite-dimensional vector space. A γ -dependent presymplectic manifold is a manifold M equipped with a smooth γ -dependent 2-form γ such that d γ = 0 for each γ ∈ . We say (M, γ ) is exact when there is a smooth γ -dependent 1-form ϑ γ such that γ = −dϑ γ .

Overview
This article contains a pair of interrelated contributions, an analog of Kruskal's nearly periodic system theory in discrete time, and an application of the new theory to structure-preserving integration of non-canonical Hamiltonian systems. Each contribution requires understanding a fair amount of technical background material to fully digest. This section therefore aims to present a (mostly) non-technical synopsis of our work. Readers interested in proceeding directly to the technical content should start with Sect. 3.
Kruskal's theory Kruskal (1962) represented an outgrowth of vigorous investigations into charged particle motion in strong magnetic fields as part of Project Matterhorn, the first serious attempt by the US government to develop nuclear fusion for peaceful energy production. Such charged particles execute rapid rotation around magnetic field lines superposed on a much slower drift motion. Moreover, the magnetic flux that threads their tight helical trajectories remains approximately constant in time over large time intervals, almost as if each particle behaved as a superconducting ring of current. Kruskal recognized that this problem comprised just one example of a rich class of multi-scale dynamical systems for which much could be said using analytical methods. These systems, now called nearly periodic, exhibit two time scales, a short one on which every trajectory is periodic, and a much longer one associated with slow drifting motion. Kruskal also recognized that the magnetic flux invariant from charged particle theory generalizes to nearly periodic systems provided the nearly periodic system admits Hamiltonian structure.
Nearly periodic systems bear similarities with the nearly integrable systems addressed by Kolmogorov-Arnol'd-Moser (KAM) theory. They depend on a small parameter that measures a degree of timescale separation. When → 0, all trajectories become very simple-even explicitly integrable in the right coordinate system. In fact, nearly periodic systems may be viewed as special examples of nearly integrable systems. However, where nearly integrable systems generically exhibit resonances nearly periodic systems never do. The limiting dynamics for a nearly integrable Hamiltonian system comprise quasiperiodic motion on invariant Lagrangian tori in phase space characterized by several fast phases. On many of these tori, the fast frequencies exhibit integer relationships that notoriously cause perturbation series to break down. But for a nearly periodic system the tori are one-dimensional, making resonance impossible. Kruskal exploited this lack of obstruction to perturbative methods to prove a remarkable result: every nearly periodic system admits a unique hidden continuous symmetry at the level of perturbation theory. While Kruskal did not prove that the series defining his symmetry converge or represent a genuine symmetry in any sense, formal existence and uniqueness proved sufficient to explain the magnetic flux invariant for charged particles, and more generally adiabatic invariants in any nearly periodic system.
A simple consequence of the existence of a formal hidden symmetry and corresponding adiabatic invariant for a nearly periodic Hamiltonian system is existence of a non-trivial class of phase space diffeomorphisms with hidden symmetries and adiabatic invariants. At a minimum, this class contains the time-t flow maps for nearly periodic Hamiltonian systems, for any value of t. The work presented in this article emerged from a desire to better understand this interesting class of maps. Following Kruskal's lead, our goal was to first define these maps axiomatically and then transcribe Kruskal's arguments to the extent possible in order to prove they have the right properties. Section 4 represents the successful culmination of this effort.
Our theory of discrete-time nearly periodic systems, or nearly periodic maps for brevity, largely parallels Kruskal's original theory in continuous time. But there is one key technical challenge that appears in discrete time and not in continuous time, which can be described as follows. Ultimately, perturbation theory for dynamical systems involves integration along unperturbed trajectories. For a continuous-time nearly periodic system, such integrations manifest themselves in the form of partial differential equations of the general type where the source function S(x, θ) and the (assumed nowhere vanishing) angular frequency ω(x) are known and f (x, θ) is the dependent variable. Provided the θ -average of S vanishes, this equation yields to solution by the method of Fourier series. No resonances appear. For nearly periodic maps, integrating along unperturbed orbits instead leads to functional equations of the form where θ 0 is a constant parameter characterizing the map's limiting rotation angle. Fourier analysis in θ reveals that solving for the n th Fourier coefficient f n requires division by 1 − exp(i n θ). This step becomes problematic whenever θ = 2 π q, where q is a rational number. Of course, this issue may be avoided formally by choosing the limiting rotation angle to be an irrational fraction of 2π . But the possibility of resonance remains in the theory, leading to a dichotomy between resonant and non-resonant nearly periodic maps. Most of Kruskal's arguments transcribe nicely to discrete time in the non-resonant case only.
In light of our recent results from Burby and Hirvijoki (2021), existence of nearly periodic maps presented interested possibilities for important applications of our theory. For example, we showed that every Hamiltonian system on an exact symplectic manifold embeds as a slow manifold for a larger nearly periodic system with a simple Hamiltonian structure. Due to the close relationship between flow maps for nearly periodic systems and nearly periodic maps, this suggests that nearly periodic maps might be useful for simulating a very broad class of Hamiltonian systems without the need for first identifying canonical variables. Section 4.2 shows that this hunch is actually correct. Given a possibly non-canonical Hamiltonian system, we first embed the dynamics as a slow manifold in a nearly periodic Hamiltonian system using the construction from Burby and Hirvijoki (2021). Then, we construct a nearly periodic map that preserves the Hamiltonian structure of the larger system and approximately integrates the nearly periodic flow. By restricting such an integrator to initial conditions that lie on the zero level set of the map's adiabatic invariant (guaranteed by the general theory of nearly periodic maps), we obtain an effective structure-preserving integrator for the original Hamiltonian dynamics. No canonical variables required.

Kruskal's Theory of Nearly Periodic Systems
In 1962, Kruskal presented an asymptotic theory (Kruskal 1962) of averaging for dynamical systems whose trajectories are all periodic to leading order. Nowadays, Kruskal's method is termed one-phase averaging (Lochak 1993), which suggests a contrast with the multi-phase averaging methods underlying, e.g., KAM theory. Since this theory provides a model for the results in this article, we review its main ingredients here. In this section only, and merely for simplicity's sake, we make the restriction = R.
Definition 2 A nearly periodic system on a manifold M is a smooth γ -dependent vector field X γ on M such that X 0 = ω 0 R 0 , where The vector field R 0 is called the limiting roto-rate, and ω 0 is the limiting angular frequency.
Remark 1 In addition to requiring that ω 0 is sign-definite, Kruskal assumed that R 0 is nowhere vanishing. However, this assumption is not essential for one-phase averaging to work. It is enough to require that ω 0 vanishes nowhere. This is an important restriction to lift since many interesting circle actions have fixed points.
Kruskal's theory applies to both Hamiltonian and non-Hamiltonian systems. In the Hamiltonian setting, it leads to stronger conclusions. A general class of Hamiltonian systems for which the theory works nicely may be defined as follows.
Definition 3 Let (M, γ ) be a manifold equipped with a smooth γ -dependent presymplectic form γ . Assume there is a smooth γ -dependent 1-form ϑ γ such that γ = −dϑ γ . A nearly periodic Hamiltonian system on (M, γ ) is a nearly periodic system X γ on M such that ι X γ γ = dH γ , for some smooth γ -dependent function H γ : M → R.
Kruskal showed that all nearly periodic systems admit an approximate U (1)symmetry that is determined to leading order by the unperturbed periodic dynamics. He named the generator of this approximate symmetry the roto-rate. In the Hamiltonian setting, he showed that both the dynamics and the Hamiltonian structure are U (1)-invariant to all orders in γ .
Definition 4 A roto-rate for a nearly periodic system X γ on a manifold M is a formal power series R γ = R 0 + γ R 1 + γ 2 R 2 + . . . with vector field coefficients such that where the second and third conditions are understood in the sense of formal power series.
Proposition 1 (Kruskal (1962)) Every nearly periodic system admits a unique rotorate R γ . The roto-rate for a nearly periodic Hamiltonian system on an exact presymplectic manifold (M, γ ) satisfies L R γ γ = 0 in the sense of formal power series.

Corollary 1
The roto-rate R γ for a nearly periodic Hamiltonian system X γ on an exact presymplectic manifold (M, γ ) with Hamiltonian H γ satisfies L R γ H γ = 0.
Proof Since [R γ , X γ ] = L R γ X γ = 0 and L R γ γ = 0, we may apply the Lie derivative L R γ to Hamilton's equation ι X γ γ = dH γ to obtain Thus, L R γ H γ is a constant function. By averaging over the U (1)-action we conclude that the constant must be zero.
To prove Proposition 1, Kruskal used a pair of technical results, each of which is interesting in its own right. The first establishes the existence of a non-unique normalizing transformation that asymptotically deforms the U (1)-action generated by R γ into the simpler U (1)-action generated by R 0 . The second is a subtle bootstrapping argument that upgrades leading-order U (1)-invariance to all-orders U (1)-invariance for integral invariants. We state these results here for future reference.
Definition 5 Let G γ = γ G 1 + γ 2 G 2 + . . . be an O(γ ) (no constant term) formal power series whose coefficients are vector fields on a manifold M. The Lie transform with generator G γ is the formal power series exp(L G γ ) whose coefficients are differential operators on the tensor algebra over M.

Definition 6 A normalizing transformation for a nearly periodic system
Proposition 2 (Kruskal) Each nearly periodic system admits a normalizing transformation.
Proposition 3 Let α γ be a smooth γ -dependent differential form on a manifold M. Suppose α γ is an absolute integral invariant for a C ∞ nearly periodic system X γ on M. If L R 0 α 0 = 0 then L R γ α γ = 0, where R γ is the roto-rate for X γ .
Proof Integral invariance means L X γ α γ = 0 for each γ ∈ . By Applying L R γ to this relationship, and using [R γ , X γ ] = 0, we obtain L X γ L R γ α γ = 0. Now let G γ be the generator of a normalizing transformation for X γ , and set Repeating this argument gives L R 0 α k = 0 for k > 1 as well. In other words L R 0 α γ = 0 to all orders in γ , which is equivalent to the theorem's claim.
According to Noether's celebrated theorem, a Hamiltonian system that admits a continuous family of symmetries also admits a corresponding conserved quantity. Therefore, one might expect that a Hamiltonian system with approximate symmetry must also have an approximate conservation law. Kruskal showed that this is indeed the case for nearly periodic Hamiltonian systems, as the following generalization of his argument shows.
Proposition 4 Let X γ be a nearly periodic Hamiltonian system on the exact presymplectic manifold (M, γ ). Let R γ be the associated roto-rate. There is a formal power series θ γ = θ 0 + γ θ 1 + . . . with coefficients in 1 (M) such that γ = −dθ γ and L R γ θ γ = 0. Moreover, the formal power series μ γ = ι R γ θ γ is a constant of motion for X γ to all orders in perturbation theory. In other words, in the sense of formal power series.
Proof To construct the U (1)-invariant primitive θ γ , we select an arbitrary primitive ϑ γ for γ and set This formal power series satisfies L R γ θ γ = 0 because whence θ γ is a primitive for γ . To establish all-orders time independence of μ γ = ι R γ θ γ , we apply Cartan's formula and Corollary 1 according to Definition 7 The formal constant of motion μ γ provided by Proposition 4 is the adiabatic invariant associated with a nearly periodic Hamiltonian system.

Nearly Periodic Maps
Nearly periodic maps are natural discrete-time analogs of nearly periodic systems.
The following provides a precise definition.
Definition 8 (nearly periodic map) Let be a vector space. A nearly periodic map on a manifold M with parameter space is a smooth mapping F : M × → M such that F γ : M → M : m → F(m, γ ) has the following properties: We say F is resonant if θ 0 is a rational multiple of 2π , otherwise F is non-resonant.
The infinitesimal generator of θ , R 0 , is the limiting roto-rate.
Differentiating this condition with respect to θ at the identity implies F * t R = R, where R denotes the infinitesimal generator for the U (1)-action. Conversely, if R is any vector field with 2π -periodic orbits and F * t R = R for each t ∈ R then the Rflow defines a U (1)-action that is a symmetry for X . Since we would like to think of nearly periodic maps as playing the part of a nearly periodic system's flow map, the latter characterization of symmetry allows us to naturally extend Kruskal's notion of roto-rate to our discrete-time setting.
Definition 9 A roto-rate for a nearly periodic map F is a formal power series . . whose coefficients are homogeneous polynomial maps from into vector fields on M such that formal power series whose coefficients are homogeneous polynomial maps from into vector fields on M. The Lie transform with generator G γ is the formal power series exp(L G γ ) whose coefficients are homogeneous polynomial maps from into differential operators on the tensor algebra over M.

Remark 2
Since the parameter γ in the previous two definitions is vector valued, the formal power series R γ and G γ may be usefully interpreted as multivariate formal power series.

Definition 11 A normalizing transformation for a nearly periodic map
Our first and most fundamental result concerning nearly periodic maps establishes the existence and uniqueness of the roto-rate in the non-resonant case. Like the corresponding result due to Kruskal in continuous time, this result holds to all orders in perturbation theory.

Theorem 1 (Existence and uniqueness of the roto-rate) Each non-resonant nearly periodic map admits a unique roto-rate.
Proof First we will show that there exists a Lie transform with generator G γ such that To that end, we introduce a convenient way of representing γ -dependent pullback operators at the level of formal power series. Let ψ γ be a smooth γ -dependent diffeomorphism on M. By the Lie derivative formula, there is a unique γ -dependent * -valued vector field W γ such that for each γ, δγ ∈ . Here ·, · denotes the natural pairing between and its dual space * . The object W γ both determines and is determined by the pullback operator ψ * γ at the level of formal power series. This follows from recursive application of the identity which may be understood as a consequence of (1) and the fundamental theorem of calculus. This can be viewed as Picard iteration of (1) or fixed point iteration of (2). The first step in the recursion is to substitute (2) with s = s 1 into (2) with s = 1, resulting in The second step involves substituting (2) with s = s 2 into the preceding expression, thereby producing a triple integral. Continuing in this manner, it is straightforward to derive the following time-ordered exponential formulas for both the pullback ψ * γ and pushforward operator ψ γ * , Upon introducing the formal power series expansion W γ = W 0 +W 1 [γ ]+W 2 [γ, γ ]+ . . . , the integrals in these formulas can be carried out, leading to the somewhat more explicit formulas The preceding discussion applies in particular to ψ * γ = F * γ . In this case, we will use the symbol V γ for W γ . The discussion also applies to the formal pullback operator as well as its inverse φ γ * = (φ * γ ) −1 . In this case we will use ξ γ in place of W γ . Thus, we have the defining identities We will now establish existence of the Lie transform with generator G γ by constructing an appropriate ξ γ . The Lie transform itself can then be constructed using the formulas (3) and (4).
If this can be done, then R γ will be a roto-rate. The equation we would like to solve is equivalent to where (φ −1 γ ) * = (φ * γ ) −1 . Formally, this is just the statement that the "diffeomorphism" γ preserves the limiting roto-rate R 0 . Instead of solving (9) directly, we will demand that its γ -derivative vanishes. This derivative condition is where V γ is readily shown to be given by Note that requiring the γ -derivative of (9) to vanish implies (9) itself since the latter is clearly satisfied when γ = 0. Also note that since F * γ is formally invertible, the derivative condition is equivalent to To solve (11), we will expand the equation in powers of γ and then argue inductively that each equation in the resulting sequence can be solved. At O(γ 0 ) we have T is any tensor field on M, and T osc = T − T , this equation is equivalent to where we have introduced the homological operator Since F γ is assumed to be non-resonant, the homological operator, regarded as a linear automorphism of the oscillating subspace of vector fields, is invertible. We may therefore solve the O(γ 0 ) equation by setting . . with polynomial degree (for ξ ) at most n − 1. Assuming the ξ k with k < n have already been determined by solving the O(γ k ) components of (11), we may therefore solve the O(γ n ) equation by setting Since we have already established that the O(γ 0 ) equation has a solution, we now conclude by induction that (11) may be solved for ξ γ to all orders in γ . It follows that a roto-rate exists.
Next we prove uniqueness of the R γ just constructed. Suppose R γ is a possibly distinct roto-rate, and consider the commutator , which implies L R 0 C 1 [γ ] = 0 by non-resonance. But since the mean of C 1 [γ ] vanishes, we must have C 1 [γ ] = 0 for all γ ∈ . The same argument applied repeatedly implies C n = 0 for all n ≥ 1. In other words, R γ and R γ commute. To finish the argument for uniqueness, we now use commutativity of R γ and R γ to find Theorem 2 (Existence of normalizing transformations) Each non-resonant nearly periodic map admits a normalizing transformation.
Proof This follows immediately from the proof of Theorem 1.

Nearly Periodic Maps With Hamiltonian Structure
As in the continuous-time theory, existence of the roto-rate leads to additional insights for nearly periodic maps that are Hamiltonian in an appropriate sense. In this subsection, we will establish the basic properties of nearly periodic maps with Hamiltonian structure. We start by defining what we mean by Hamiltonian structure.

Definition 12 (Presymplectic nearly periodic map) A Presymplectic nearly periodic
Definition 13 (Hamiltonian nearly periodic map) A Hamiltonian nearly periodic map on a γ -dependent presymplectic manifold (M, γ ) is a nearly periodic map F such that there is a smooth (t, γ )-dependent vector field Y t,γ with the following properties: Lemma 1 Each Hamiltonian nearly periodic map is a presymplectic nearly periodic map.
the last identity may be rewritten as the sequence of identities e ikθ 0 k 0 = k 0 , k ∈ Z. But by non-resonance of F, 1 − e ikθ 0 is nonvanishing for each k. We conclude that k 0 = 0 for nonzero k, or, equivalently, L R 0 0 = 0. Presymplecticity of F for nonzero γ implies F * γ γ = γ for each γ ∈ . Applying the Lie derivative L R γ to this identity and using F * In other words, α γ = L R γ γ is (formally) an absolute integral invariant for F γ . By the argument from the previous paragraph, we see immediately that α 0 = 0. To finish the proof, we will use integral invariance together with existence of a normalizing transformation to find that α γ = 0 to all orders in γ . This argument will parallel the proof of Proposition 3.
Using presymplecticity of the roto-rate, we may now use a version of Noether's theorem to establish existence of adiabatic invariants for many interesting presymplectic nearly periodic maps.
Theorem 4 (Existence of an adiabatic invariant) Let F be a non-resonant presymplectic nearly periodic map on the exact γ -dependent presymplectic manifold (M, γ ) with roto-rate R γ . Assume one of the following conditions is satisfied.
(2) M is connected and the limiting roto-rate R 0 has at least one zero.
There exists a smooth γ -dependent 1-form θ γ such that L R γ θ γ = 0 and −dθ γ = γ in the sense of formal power series. Moreover the quantity satisfies F * γ μ γ = μ γ in the sense of formal power series. In other words, μ γ is an adiabatic invariant for F. Proof By Theorem 3, a primitive θ γ with the desired properties may be constructed as in the proof of Proposition 4.
To establish adiabatic invariance of μ γ , first we compute the exterior derivative of μ γ using Cartan's formula to obtain dμ γ = ι R γ γ . Since both R γ and γ are F γ -invariant, it follows that dF * γ μ γ = dμ γ . The difference c γ ≡ F * γ μ γ − μ γ must therefore be a formal power series whose coefficients are homogeneous polynomial maps from into locally constant functions on M. To complete the proof, we must demonstrate that c γ = 0 to all orders.
First suppose F is Hamiltonian. Then, there is a smooth (t, γ )-dependent Hamiltonian vector field Y t,γ with Hamiltonian H t,γ whose t = 1 flow is equal to F γ . Let F γ t denote the time-t flow map for Y t,γ with F γ 0 = id M . By the fundamental theorem of calculus, the definition of Lie derivative, and Cartan's formula, we therefore have Applying exp(θ L R 0 ) to this identity and averaging in θ gives the desired result. Finally suppose that M is connected and that R 0 (m 0 ) = 0 for some m 0 ∈ M. Let G γ be the generator of a normalizing transformation. We have It follows that c γ is zero on the connected component of M containing m 0 . But since M is connected, c γ is therefore zero everywhere, as claimed.

Remark 3 A simple example illustrates how existence of an adiabatic invariant can fail.
Let M = S 1 × R (ζ, I ), γ = dζ ∧ d I = −d(I dζ ), and = R. The mapping F(ζ, I , γ ) = (ζ + θ 0 , I + γ ) defines a non-resonant nearly periodic map for almost all θ 0 ∈ U (1). The roto-rate is given to all orders by R γ = ∂ ζ . Moreover, F γ is areapreserving for each γ , and hence presymplectic. The quantity μ γ = ι R γ (I dζ ) = I from (14) is clearly not an adiabatic invariant for F since F * γ I = I + γ . Note that, in this example, the R 0 -orbits are not contractible and that F is presymplectic but not Hamiltonian.

Geometric Integration of Noncanonical Hamiltonian Systems Using Nearly Periodic Maps
Let defines a metric tensor on Q. Equip the tangent bundle π : T Q → Q with the "magnetic" symplectic form * = π * ω + , where is a real parameter and is the pullback of the canonical symplectic form on T * Q along the bundle map T Q → T * Q defined by g. We may also define a natural Hamilton function on T Q, As explained in detail in Burby and Hirvijoki (2021), H * defines a Hamiltonian nearly periodic system whose slow manifold dynamics recovers the dynamics of H on Q as → 0. The limiting roto-rate is where X H denotes the Hamiltonian vector field on Q associated with H , and the angular frequency function ω 0 = 1. Moreover, the adiabatic invariant associated with H * ensures that the slow manifold enjoys long-term normal stability. It is crucial that the metric g is determined by an almost complex structure J compatible with ω for these results to hold. If g is a more general metric tensor then the larger system on T Q need not be nearly periodic, and an adiabatic invariant need not exist. The purpose of this section is to combine observations from Burby and Hirvijoki (2021) with the theory of nearly periodic maps in order to construct a geometric numerical integrator for H . The integrator will be given as an implicitly defined mapping on T Q that is provably presymplectic nearly periodic with limiting roto-rate R 0 . We will show that this mapping admits a slow manifold diffeomorphic to Q on which iterations of the map approximate the H -flow. In fact, the mapping is a slow manifold integrator in the sense described in Burby and Klotz (2020). In addition, we will argue using the Noether theorem for nearly periodic maps that this discrete-time slow manifold enjoys long-term normal stability. This ensures that the mapping on T Q will function effectively and reliably as a structure-preserving integrator for the original Hamiltonian system on Q. We remark that the results described in this section provide a general solution to the problem of structure-preserving integration of non-canonical Hamiltonian systems on exact symplectic manifolds. For a completely different approach that is less geometric, we refer readers to Kraus (2017).
We begin with some preliminary remarks.
Remark 4 It will be convenient to work with the parameter = √ instead of . There are technical reasons for doing so that will not be discussed here; however, an obvious physical benefit will be that may be interpreted as a timestep. The symplectic form on T Q will therefore be given by Remark 5 It is useful to describe the goal of our construction in more concrete terms. We aim to find a smooth -dependent mapping : T Q → T Q that is both nonresonant nearly periodic with limiting roto-rate R 0 given by (15) and symplectic, * * = * , for all 1. Since Q is connected and R 0 has a manifold of fixed points of the form {(q, X H (q))} ⊂ T Q, Theorem 4 (Noether's theorem for nearly periodic maps) ensures that the mapping we seek will admit an adiabatic invariant μ .

Remark 6
We may determine the leading-order term in the formal power series μ = μ 0 + μ 1 + 2 μ 2 +. . . using only the explicit expressions for * and R 0 in conjunction with the general existence theorem (Theorem 4). Recall that the theorem says that the adiabatic invariant is given by μ = ι R , where R is the roto-rate and is a U (1)-invariant primitive for * . In particular, the roto-rate must satisfy Hamilton's equation ι R * = dμ with Hamiltonian μ . We apply this theorem as follows.
We have now established that the adiabatic invariant for the nearly periodic map we aim to construct must have the form μ = 4 μ 4 + 5 μ 5 + . . . . We can determine an explicit expression for μ 4 as follows. By the above remarks, we must have R = (J ∂ v μ 4 ) · ∂ v + l.o.t., which implies in particular that R 0 = (J ∂ v μ 4 ) · ∂ v . Since the desired form of R 0 is known, we therefore obtain the following partial differential equation for μ 4 :

The general solution of this equation is given by
is an arbitrary function of q. To determine χ , we evaluate the formula μ = ι R at a fixed point (q, v (q)) to find 0 = lim →0 −4 μ (q, v (q)) = μ 4 (q, v 0 (q)) = χ(q). Note that we have used the formula for v 0 (q) = X H (q) for fixed points of R 0 . We conclude that the adiabatic invariant must have the general form This formula will be useful later when we argue for long-term normal stability.
To construct the mapping : T Q → T Q, we begin by introducing the mixedvariable generating function where (q, q) = (q/2 + q/2, q − q) is a diffeomorphism Q × Q → T Q, and : T Q → R is given by Here, ·, · is shorthand for g(·, ·), the integral is taken along the straight line connecting q with q, X H = −J∇ H is the Hamiltonian vector field associated with H , and θ 0 ∈ /{0, π}. The metric tensor, the Hamiltonian, and the Hamiltonian vector field are evaluated at the midpoint x = (q + q)/2. The variables q and q should be interpreted as the "old" and "new" points in the symplectic manifold Q. While it is not necessary to go into the details here, the function S may be understood as an approximation of Jacobi's solution of the Hamilton Jacobi equation for the Hamiltonian H * .

Definition 14
The symplectic Lorentz map is the mapping : T Q → T Q : (q, v) → (q, v) defined by the implicit relations Here g q (v, dq) is the linear map T q Q → R given by w q → g q (v, w q ), and g q (v, dq) is defined analogously. (q, v, ) for in a neighborhood of 0 ∈ R. Moreover, it preserves the -dependent symplectic form * and satisfies

Proposition 5 The symplectic Lorentz map is well-defined and smooth in
Proof First we will construct a convenient moving frame on Q × Q onto which we will resolve the implicit relations (18) and (19). We will start by building a frame on T Q and then finish by pulling back along the mapping : Q × Q → T Q defined above. Without loss of generality, assume Q = R n for an even integer n and let (x i , ξ i ) denote the standard linear coordinate system on T Q. Fix (x, ξ) ∈ T Q and let γ : [−1, 1] → Q be a smooth curve in Q with γ (0) = x. Relative to the Riemannian structure defined by the metric g there is a unique horizontal lift γ : [−1, 1] → T Q withγ (0) = (x, ξ). In this manner, to each (x, ξ) ∈ T Q and each tangent vector w ∈ T x Q we assign a lifted tangent vectorw ∈ T (x,ξ ) T Q. Applying this construction point-wise to the coordinate vector fields ∂ x i on Q, we obtain linearly independent vector fields∂ x i on T Q. The collection of 2n vector fields (∂ x i , ∂ ξ i ) comprise a frame on T Q. A frame on Q × Q is then furnished by the vector fields the vector fields∂ x i may be written as Therefore, we may write the following explicit formulas for the frame (U i , A i ): where we remind the reader that x = (q + q)/2 and ξ = q − q. Next, we rewrite (18) and (19) as a single equation on Q × Q, and then take components along the frame (U i , A i ) to obtain Here, we have applied the useful formulas To show that these implicit equations define a smooth -dependent mapping , we first introduce the new variable = −2 (ξ − X H (x)) and then observe that when expressed in terms of the implicit equations above may be written in the form Note in particular that the term 1 12 2 ∂ k ω ji X k H X j H exactly cancels the second-order part of 1 0 λ − 1 2 ξ j ω ji (λ) dλ. Dividing these expressions by 2 implies that there are smooth functions Z 1 , Z 2 : R n × R n × R n × R n × R → R n given by such that the implicit equations defining the symplectic Lorentz map are satisfied if and only if When = 0, the unique solution of these equations for and v is which is invertible. The implicit function theorem therefore implies there is a unique pair of smooth functions (x, v, ), v(x, v, ) defined in an open neighborhood of Since is related to q by another simple application of the implicit function theorem establishes existence and smoothness of the symplectic Lorentz map : (q, v) → (q, v). We have also shown that 0 has the desired form (q, v) → (q, v 0 ). Symplecticity of now follows immediately from applying the exterior derivative to (21).

Corollary 2 The symplectic Lorentz map is a presymplectic nearly periodic map. It is non-resonant provided
We have thus constructed a nearly periodic map with the desired roto-rate and integral invariant. Now we must establish the precise sense in which the symplectic Lorentz map , which is a priori a mapping T Q → T Q, functions as a consistent numerical integrator for the Hamiltonian system X H on Q. The first hint as to how this might work is that the limit map 0 admits a manifold of fixed points given by This limiting invariant manifold, being the graph of X H , is manifestly diffeomorphic to Q. Thus, if 0 can be continued to an invariant manifold for with = 0, we would automatically obtain dynamics on Q that could be compared with those of X H by restricting to . Unfortunately, 0 is unlikely to continue as a true invariant manifold since each fixed point on 0 is of elliptic type. Instead, we can obtain the following weaker result. Roughly speaking, it says that there is a unique invariant continuation of 0 at the level of formal power series in .
Proof Expanding the condition (22) in powers of leads to an infinite sequence of constraints that the formal power series v * must obey. Simultaneous satisfaction of each constraint in the sequence is equivalent to (22). The first two constraints are given explicitly by where we have introduced the formal series expansions More generally, the n th equation in the sequence has the form where s n (q) depends only on coefficients of the power series expansion for and coefficients v * k with k < n. Invertibility of L(q) therefore implies that there is a unique formula for v * n for each n. The formal power series v * defined in this manner satisfies (22) by construction.
So while we do not obtain a genuine invariant manifold diffeomorphic to Q, we do obtain a family of approximate invariant manifolds diffeomorphic to Q given by truncations of the formal power series v * . Using arguments comparable to those presented in Burby and Hirvijoki (2021), it is possible to show that truncations of v * may be constructed so their graphs agree with the zero level set of the adiabatic invariant μ for to any desired order in . Adiabatic invariance of μ can then be used to prove the existence of manifolds (n) close to 0 with the following schematic normal stability property: • For each N > 0 and (q, v) within α(n) of (n) the point k (q, v) remains within β(n) of (n) for k = O( −N ). Here α and β are monotone increasing functions of n.
We will not attempt to prove such a result in full generality here. However, since we have already determined the form of the leading term in the adiabatic invariant series (in this case μ = 4 μ 4 + O( 5 )), we can prove a special case of the general result without much effort. First we establish the timescale over which μ 4 is well-conserved.
Definition 15 Given a compact set C ⊂ T Q, a point (q, v) is positively contained if k (q, v) ∈ C for all nonnegative integers k and all in a neighborhood of 0.

Proposition 7 For each N > 0 and compact set C ⊂ T Q there is a positive,independent constant M such that
Proof First we will obtain a useful estimate for the degree of conservation of an arbitrary truncation of the adiabatic invariant series. Let μ = −4 μ = μ 4 + μ 5 + . . . denote the reduced adiabatic invariant for the symplectic Lorentz map. Define ) and μ is -invariant to all orders in , there is a constant M N , depending on both C and N , such that For positively contained (q, v), we may apply this formula repeatedly to obtain an estimate for the change in μ (N ) after k positive timesteps, Next, we draw implications from the previous inequality together with a bound on the difference between μ (0) = μ 4 and μ (N ) . There must be another positive constant M N , depending on both C and N , such that In light of the inequality (23), this implies that for each positively contained (q, v) the change in μ (0) after k positive timesteps is at most Apparently, the change in μ 4 = μ (0) is at most O( ) as long as k N +1 = O( ). We therefore obtain the desired inequality with k * ( , N ) = −N − 1.

Remark 7
If "positively contained" fails, then the next best thing would be having uniform bounds (in T Q and ) on the derivatives of . The precise form of these bounds would depend on the details of the underlying Hamiltonian system on Q. Then, the proof of the previous proposition would go through practically unchanged. In the absence of both positively containment and uniform boundedness, things get trickier, and we have no general answer as to the validity of Proposition 7.
Using this result together with the explicit form of μ 4 , we now easily obtain the following normal stability result for the almost invariant set given by the graph of X H .

Proposition 8 Let C ⊂ T Q be a compact set and set
Let | · | denote the velocity norm provided by the metric tensor g. For each N > 0, V 0 > 0, and positively contained (q, v) ∈ C that satisfies there is a positive constant V 1 such that Proof Let (q, v) ∈ C be positively contained and suppose |v − X H (q)| q < V 0 . By Proposition 7, we have for some N -dependent constant M and k ∈ [0, k * ( , N )]. But since μ 4 (q, v) = 1 2 |v − X H (q)| 2 q we can apply this inequality to obtain Taking a square root gives the desired result.
In the above sense, the graph of X H behaves much like a true invariant set over very large time intervals. Of course, the invariance need not be exact, but may include oscillations around the graph of amplitude √ . The amplitude of these oscillations can be reduced by considering manifolds that better approximate the zero level set of μ , but, as mentioned earlier, we will not pursue this matter further in this article.
To complete the picture of how the symplectic Lorentz map may be used as an integrator for X H on Q, we will now describe the precise sense in which the map's dynamics approximate the H -flow. We start with a simple estimate that says the qcomponent of the symplectic Lorentz map approximates the time-flow of X H on Q with an O( 5/2 ) error, provided the map is applied in an O( 1/2 ) neighborhood of the graph {v = X H (q)}.
Proof In the proof of Proposition 5, we already established and where x = q/2+q/2. Implicit differentiation of these formulas together with Taylor's theorem with remainder therefore implies Combining this result with our earlier estimate of the normal stability timescale for {v = X H (q)} in Proposition 8 finally allows us to conclude that the q-component of the symplectic Lorentz map provides a persistent approximation of the H -flow over very large time intervals provided initial conditions are chosen close enough to the graph {v = X H (q)}.
Corollary 3 (Persistent approximation property) Let C be a compact set and let (q, v ) ∈ C be a smooth -dependent point in C that is positively contained for Proof Proposition 8 ensures that the iterates (q k , v k ) remain within O( 1/2 ) of {v = X H (q)} for k in the desired range. Thus, Proposition 9 applies to each iterate individually, which is precisely the desired result.
In summary, we have established the following remarkable properties of the symplectic Lorentz map .
1. It is symplectic on T Q, when T Q is endowed with the magnetic symplectic form * = π * ω + 2 .
2. Its q-component provides an approximation of the time-flow of X H with O( 5/2 ) local truncation error when applied to points in T Q within O( 1/2 ) of {v = X H (q)}. 3. If an initial condition is chosen to lie within 1/2 of {v = X H (q)} then it will remain within 1/2 of the same set for a number of iterations that scales like −N for any N .

Hidden-Variable Newtonian Gravity
In this section, we will use nearly periodic maps to construct a discrete-time model of Newtonian gravitation where the gravitational constant has a dynamical origin. Let defines a continuous-time nearly periodic Hamiltonian system with equations of motionṗ The angular frequency function is ω 0 = 1, the limiting roto-rate is R 0 = −q ∂ p + p ∂ q , and the corresponding U (1)-action is θ (q, p, Q, P) = (cos θ q + sin θ p, cos θ p − sin θ q, Q, P). When = 0, the system's flow is F t (q, p, Q, P) = t (q, p, Q, P). Intuitively, the (q, p) variables correspond to a fast oscillator that couples nonlinearly to a mechanical system parameterized by ( Q, P). The averaged Hamiltonian for the coupled system is where μ 0 = 1 2 ( p 2 + q 2 ) is the leading-order adiabatic invariant. We therefore expect the slow variables ( Q, P) to behave like a particle in d-dimensional space subject to the effective potential V ( Q) + μ 0 W ( Q).
We will construct a Hamiltonian non-resonant nearly periodic map that accurately simulates the slow dynamics for this system while "stepping over" the shortest scale 2π/ω 0 ∼ 1. If h ∈ R denotes the temporal step size, these requirements translate into symbols as 1 h −1 . Upon introducing the parameters δ = 1/h, = h, and γ = ( , δ), we may state our requirement equivalently as |γ | 1. Our construction will now proceed using the method of mixed-variable generating functions.
The exact Type I generating function for this problem can be characterized by Jacobi's solution of the Hamilton-Jacobi equation, which is given by where (q(t), Q(t), p(t), P(t)) satisfies Hamilton's equations, and the boundary conditions q(0) = q, Q(0) = Q, q(h) = q, Q(h) = Q. In the setting of variational integrators (Marsden and West 2001), this is referred to as the exact discrete Lagrangian, and there are also exact discrete Hamiltonians (Leok and Zhang 2011) corresponding to Type II and Type III generating functions. One possible way to construct a computable approximation of the exact Type I generating function is to observe that it can also be expressed as Then, one can construct a computable approximation by replacing the infinitedimensional function space C 2 ([0, h], M) with a finite-dimensional subspace, and replacing the integral with a numerical quadrature formula, which yields a Galerkin discrete Lagrangian. Under a number of technical assumptions, the resulting variational integrators -converge to the exact flow map (Müller and Ortiz 2004), and a quasi-optimality result (Hall and Leok 2015) implies that the rate of convergence is related to the best approximation properties of the finite-dimensional function space used to construct the Galerkin discrete Lagrangian. In general, this means that a good integrator can be constructed by choosing a finite-dimensional function space that is rich enough to approximate the exact solutions well, and using a quadrature rule that is accurate for that choice of function space. This might entail augmenting the function space with the solution of the fast dynamics when the slow variables are frozen, and then using a quadrature rule that is well adapted to highly oscillatory integrals, like Filon quadrature (Iserles and Nørsett 2004). In this case, the problem exhibits a fast-slow structure that lends itself to a hybrid approximation. We exploit the timescale separation to approximate the fast variables of the dynamics (q(t), p(t)) by the exact solution of the = 0 limiting system, where the slow variables ( Q(t), P(t)) are frozen, which leads to a sinusoidal solution for (q, p). Furthermore, because the timestep h is assumed to be large enough that the fast variables perform many revolutions in that time, we anti-alias the fast dynamics by replacing the revolutions by just the fractional part of the revolutions, which we denote by θ 0 , and which is assumed to be some irrational multiple of 2π , so that the invariant distribution remains the same. The component of the action integral associated with the fast variables can be evaluated analytically in this case. As for the slow variables, we adopt an approach that can be used to derive the implicit midpoint rule, which is a symplectic integrator for Hamiltonian systems. This involves approximating the solution space by linear functions, so Q(t) is uniquely determined by the boundary conditions, and approximating the integral by the midpoint rule. The use of mixed quadrature approximations of the action integral was the basis for implicit-explicit variational integrators for fast-slow systems (Stern and Grinspun 2009).
Note that ω 0 = 1, then for the (q, p) dynamics to have a θ 0 rotation in time h, the solution is given by where θ 0 /h = 1. Since q(t) = sin(t) p + cos(t)q, then q = q(h) = sin(θ 0 ) p + cos(θ 0 )q, and hence, p = q−cos(θ 0 )q sin(θ 0 ) . Therefore, the (q, p) dynamics, expressed in terms of the boundary data, is given by For the slow degrees of freedom, we consider a linear interpolant in time, Q(t) = Q+ Q− Q h t, from which the momentum becomes By using = h, replacing q h 2 with q(0) = q, and combining it with the first term coming from the fast dynamics, we obtain the following Type I generating function: where θ 0 is some irrational multiple of 2π . The implicit relations define a γ -dependent symplectic map F γ : (q, p, Q, P) → (q, p, Q, P). For small γ , we claim this map accurately captures the averaged dynamics of the slow variables in the system (25)-(28) and preserves the adiabatic invariant μ 0 over very large time intervals. To show this, we first compute the derivatives in (33) to explicitly write the defining equations for F γ as The first pair of equations can be solved explicitly for p and q, giving Adding and subtracting the last P and P equations then gives These formulas show that F 0 = θ 0 , which implies that F γ comprises a non-resonant, Hamiltonian, nearly periodic map. In particular, this map admits an all-orders adiabatic Moreover, Eqs. (36) and (37) provide a consistent numerical scheme for the averaged dynamics of the slow variables. To see this, observe that the average of q 2 in (36) after many iterations tends to μ 0 , which implies that, on average, (36) and (37) comprise the implicit midpoint scheme applied to the continuous system's averaged dynamics. Note that the relationship between the physical timestep h and is h = / . Also note that the approximation q(h/2) ≈ q(0) used when approximating the action integral is not systematic due to the rapid oscillations in q(t). This was done merely for the sake of obtaining an especially simple time advance. A more systematic approach would adopt Filon-type quadrature for the part of the integrand involving both slowly and rapidly varying terms, but the resulting nearly periodic map would have the same qualitative properties as the one introduced here. A planar N -body problem in Cartesian (x, y)-coordinates provides a convenient sandbox for testing the novel scheme (34)-(37). Assume two bodies, labeled by the position vectors Q 1 = (Q 1,x , Q 1,y ) and Q 2 = (Q 2,x , Q 2,y ) and the respective momentum vectors P 1 = (P 1,x , P 1,y ) and P 2 = (P 2,x , P 2,y ), to orbit an infinitely massive body at the origin. The potential V ( Q) is therefore Also assume the two bodies to interact via the additional central potential The instantaneous value of q 2 therefore indicates the strength of the coupling of the two bodies via the temporal evolution of the -perturbed (q, p) oscillator. The behavior of the scheme (34)-(37) is illustrated in Fig. 1 together with the numerical solution from the well-known implicit midpoint scheme which is symplectic for canonical Hamiltonian systems and generally considered a good scheme for stiff problems. For both integrators, we set the system parameter to = 0.001. In Fig. 1, the columns (a), (b), and (c) correspond to implicit midpoint scheme with time steps h = 0.1, h = 4.0, and h = 100.0, respectively. The columns (d) and (e) correspond to the fast-slow scheme with a time step of h = 100.0 and the angle variable being (d) non-resonant θ 0 = 2.0 and (e) resonant θ 0 = π . The column (a) can be considered as the reference solution, which the non-resonant fast-slow integrator in column (d) closely matches.
although care is needed in choosing the saturation value for phase angle for the limit of the nearly periodic map.

Reduced Guiding-Center Motion
We now apply the general theory developed in Sect. 4.2 to motion of a charged particle in a strong magnetic field of the special form B(x, y, z) = B(x, y) e z , where (x, y, z) denotes the usual Cartesian coordinates on R 3 and B is a positive function. Let q = (x, y) ∈ Q = R 2 and introduce a symplectic form ω on Q using the formula Here, the components of the 1-form α may be interpreted as the physicist's vector potential for B. Also define the Hamiltonian function H : Q → R according to where μ is a positive constant parameter. The corresponding Hamiltonian vector field is given by where R π/2 is the rotation matrix R θ evaluated at π/2, R θ = cos θ − sin θ sin θ cos θ .
Physically, this Hamiltonian vector field describes the motion of a charged particle's guiding center (Northrop 1963) (x, y). The parameter μ is the magnetic moment, and X H is also known as the ∇ B-drift velocity. We remark that readers familiar with the Hamiltonian formulation of guiding center motion Cary and Brizard 2009) may be used to seeing these equations derived from the Lagrangian L : T Q → R given by We also remark that in this formulation of guiding center dynamics we have used translation invariance along z to eliminate the (constant) velocity along the magnetic field and the corresponding ignorable coordinate z.
In order to construct the symplectic Lorentz map for this system, we begin by observing that the complex structure which can be iterated for η and ξ . After that, one solves for v by evaluating, for example, the expression Next, we perform some numerical tests. First, we choose a magnetic field where α introduces a small perturbation to the otherwise constant magnetic field. For the original systemq = X H (q), this field results in circular orbits for q. We then investigate the solutions of the symplectic Lorentz map and compare them with the classic RK4 integrator and the implicit midpoint applied to the original system. For the scaling of time, we choose τ = α −1 . Choosing an initial point q = (1, 1), parameters B 0 = 1, μ = 1.0, α = 0.001, and initializing the Lorentz map with v = X H (q), we run the simulation for 60'000 steps of size = 0.1. This is enough to demonstrate the Fig. 2 Comparison of the guiding-center RK4 integrator, implicit midpoint, and the symplectic Lorentz map in the simple magnetic field case. The orbit radius |q| of the RK4 integrator deteriorates clearly while the symplectic Lorentz map and the implicit midpoint manage to retain the radius within stable limits deterioration of the RK4 method while the symplectic Lorentz map and the implicit midpoint preserve the orbit in place, as seen in Fig. 2. The average number of iterations for solving the discrete system of equations and the average total execution times for the different algorithms are recorded in Table 1. Solving the nonlinear equations to machine precision results on average 16 iterations for the implicit midpoint and 22 iterations for the symplectic Lorentz map during every time step. The execution time of the Lorentz map is approximately six times that of the implicit midpoint method, credited to the larger number of iterations required and the need for additional quantities to be evaluated, such as the line integrals present in the generating function. This limited comparison with other methods should not be construed as the last word on the subject. The important takeaways include (a) the number of implicit iterations needed for one step of the symplectic Lorentz map is comparable to the number of iterations required for implicit midpoint applied to the original problem, even though the number of function evaluations is higher for symplectic Lorentz, (b) the symplectic Lorentz map achieves similar long-time solution quality as a popular scheme for integrating non-dissipative systems, and (c) the symplectic Lorentz map has provable structure-preserving properties, while implicit midpoint does not in this case. (Implicit midpoint is known to be symplectic for canonical Hamiltonian systems, but not in general for non-canonical systems like the one considered here.) Next, we consider the magnetic field  (6000) have been chosen to illustrate the non-trivial beating structure of the adiabatic invariant whose level sets have a "figure-eight" structure. By energy conservation, the guidingcenter orbits should reflect this pattern. Choosing a time step of = 0.05 and τ = 1.0, we run the symplectic Lorentz map for 6'000 steps and illustrate both the orbits and the evolution of the postulated adiabatic invariant Fig. 3. The orbits appear stable and well-confined to their respective phase-space domains, and the adiabatic invariant remains within bounds while oscillating with a non-trivial beating structure.
For the same magnetic field, we performed a pair of tests that probe the robustness of the discrete-time adiabatic invariant μ. We introduce an empirical estimate of the breakdown time for μ conservation and compute how that estimate varies with the parameters θ 0 and . Our estimate is based on the observation that μ typically oscillates about a time-varying mean value μ with an approximately constant oscillation amplitudeμ. We estimate that breakdown has occurred after n iterations when μ(n )−μ(0) >μ. We then define the breakdown time estimate to be T breakdown = n . Results from our sensitivity studies are displayed in Fig. 4. While the general theory predicts that the breakdown time should scale as fast as −N for any nonnegative integer N , the observable asymptote in T breakdown ( ) appears well-approximated by −3.5 . We presently lack understanding of the origin of the scaling exponent −3.5. The theory also predicts that adiabatic invariance should be less robust when θ 0 /2π is rational. This prediction is consistent with the plot of T breakdown (θ 0 ), which shows intermittent depressions in the breakdown time superposed on a strong upward trend as θ 0 approaches π . We hypothesize that these depressions occur at small denominator rational values of θ 0 /2π that produce nonlinear self-resonance in the integrator. As with the scaling exponent, we presently lack detailed understanding for the dramatic increase in observed breakdown time as θ 0 approaches π . Colorscale indicates value of θ 0 , ranging from θ 0 = π/4 (purple) to θ 0 = 3π/4 (red). Initial condition is x = 2.0, y = 0.0. Central dark green line is −3.5 , for reference. Theory predicts that the breakdown time should scale like −N for any positive N when is small enough. While superpolynomial scaling of the breakdown time as a function of is not apparent in these computations, it cannot be ruled out given the limited range of values considered. Computing the breakdown time for appreciably smaller values of rapidly becomes prohibitively expensive because the adiabatic invariant is so well conserved

Discussion
In this article, we have introduced and developed the theoretical foundations of nearly periodic maps. These maps provide a discrete-time analog of Kruskal (1962) continuous-time nearly periodic systems. The limiting dynamics of both nearly periodic systems and nearly periodic maps translate points along the orbits of a principal circle bundle. In the continuous-time case, each limiting trajectory ergodically samples an orbit. In discrete time, non-resonance appears as an additional requirement for ergodic sampling. As a first major application of nearly periodic maps, we used them to construct a class of geometric integrators for Hamiltonian systems on arbitrary exact symplectic manifolds.
Kruskal's principal interest in continuous-time nearly periodic systems came from their relationship to the theory of adiabatic invariants. In the paper (Kruskal 1962), Kruskal showed that nearly periodic systems necessarily admit approximate U (1)symmetries. He then went on to deduce that this approximate symmetry implies the existence of an adiabatic invariant when the underlying nearly periodic system happens to be Hamiltonian. The theory of nearly periodic maps is satisfying in this respect since it establishes the existence of a discrete-time adiabatic invariant for nearly periodic maps with an appropriate Hamiltonian structure. Moreover, the arguments used in the existence proof parallel those originally used by Kruskal. (See Thm. 4.) It is useful to place the integrators developed in this article in the context of previous attempts at geometric integration of non-canonical Hamiltonian systems. Based on the observation (Arnold 1989) that Hamiltonian systems on exact symplectic manifolds admit degenerate "phase space Lagrangians" (Cary and Littlejohn 1983), Qin and Guan (2008) proposed direct application of the theory of variational integration (Marsden and West 2001) to phase space Lagrangians for non-canonical systems. While initial results looked promising, further investigations by Leland Ellison (2016); Ellison et al. (2017) revealed that the most intuitive variational discretizations of phase space Lagrangians typically suffer from unphysical instabilities known as "parasitic modes" (Hairer et al. 2006). As noticed first in Rowley and Marsden (2002), the origin of these parasitic modes is related to a mismatch between the differing levels of degeneracy in the phase space Lagrangian and its discretization. Our integrators may be understood as modifications of those studied by Qin and Ellison that stabilize the parasitic modes over very large time intervals by way of a discrete-time adiabatic invariant. This "adiabatic stabilization" mechanism is conceptually interesting since it suppresses numerical instabilities without resorting to the addition of artificial dissipation. Also of note, adiabatic stabilization differs from the stabilization mechanism proposed by Ellison in Ellison et al. (2017), wherein the phase space Lagrangian is discretized so that it has the same level of degeneracy as its continuous-time counterpart. While Ellison's "properly degenerate" discretizations apply to a very limited class of noncanonical Hamiltonian systems, (see Ellison et al. 2017 for the precise limitations) the adiabatic stabilization method discussed here applies to any Hamiltonian system on an exact symplectic manifold.
In the preprint (Kraus 2017), Kraus has developed an alternative approach to structure-preserving integration of non-canonical Hamiltonian systems based on projection methods. In contrast to our approach, this technique is designed to produce integrators that preserve the original system's symplectic form, rather than a symplectic form on a larger space. However, there is no geometric picture for why Kraus' method ought to have this property. In fact, Kraus finds that geometrically reasonable variants of his method are not symplectic. The structure-preserving properties of our method are easier to understand in this respect, since they follow from the standard theory of mixed-variable generating functions for symplectic maps. Both techniques warrant further investigation.
As a final remark concerning relationships between the theory developed here and previous work, it is worthwhile highlighting the technique introduced by Tao in Tao (2016) for constructing explicit symplectic integrators for non-separable Hamiltonians. The latter technique applies to canonical Hamiltonian systems with general Hamiltonian H (q, p). It proceeds by constructing a canonical Hamiltonian system in a space with double the dimension of the original (q, p) space, and then applying splitting methods to the larger system. Much like the symplectic Lorentz system introduced in Burby and Hirvijoki (2021), and exploited in Sect. 4.2, Tao's larger system contains a copy of the original system as a normally elliptic invariant manifold. This suggests that Tao's construction might be interpreted as an application of nearly periodic maps. It is a curious fact, however, that Tao's error analysis suggests the oscillation frequency around the invariant manifold cannot be made to be arbitrarily large. This indicates nearly periodic map theory is not an appropriate tool for understanding Tao's results. It would be interesting to investigate whether or not nearly periodic map theory can be used to sharpen Tao's estimates.